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Abstract 

Multiple  sensor  arrays  distributed  over  a  region  provide  the  means  for  accurate  localization  of 
the  {x^y)  position  of  a  source.  When  microphone  arrays  are  used  to  measure  aeroacoustic  signals 
from  ground  vehicles,  random  fluctuations  in  the  air  lead  to  frequency-selective  coherence  of  the 
signals  that  arrive  at  widely-separated  arrays.  We  have  shown  previously  that  even  in  cases  of 
imperfect  spatial  coherence,  improvements  in  source  localization  accuracy  are  possible  when  the 
data  from  widely-separated  arrays  are  processed  jointly  by  a  fusion  center.  Further,  we  have 
shown  that  a  distributed  processing  scheme  involving  bearing  estimation  at  individual  arrays 
and  time-delay  estimation  between  pairs  of  widely-separated  sensors  performs  nearly  as  well  as 
the  optimum  scheme,  with  significantly  lower  communication  bandwidth.  These  results  were 
obtained  by  studying  the  Cramer-Rao  bound  (CRB)  on  source  localization  accuracy  based  on  a 
statistical  model  for  the  data  measured  at  the  sensors.  The  contributions  of  this  paper  include 
the  presentation  of  more  accurate  performance  bounds  (Ziv-Zakai),  and  the  development  of  a 
narrowband  subspace  algorithm  for  source  localization  with  distributed  arrays  and  partially 
coherent  signals.  We  demonstrate  through  analysis,  simulation,  and  processing  of  measured 
data  that  the  performance  of  both  algorithms  is  limited  by  ambiguities  that  arise  from  the 
narrowband  signals  and  the  large  spacing  between  arrays. 


1  Introduction 

We  are  concerned  with  estimating  the  location  (x,,  j/^)  of  a  wideband  source  using  multiple  sensor 
arrays  that  are  distributed  over  an  area.  We  consider  schemes  that  distribute  the  processing  between 
the  individual  arrays  and  a  fusion  center  in  order  to  limit  the  communication  bandwidth  between 
arrays  and  fusion  center.  Triangulation  is  a  standard  approach  for  source  localization  with  multiple 
sensor  arrays.  Each  array  estimates  a  bearing  and  transmits  the  bearing  to  the  fusion  center,  which 
combines  the  bearings  to  estimate  the  source  location  (a:^,  y»).  Triangulation  is  characterized  by  low 
communication  bandwidth  and  low  complexity,  but  it  ignores  coherence  that  may  be  present  in  the 
wavefronts  that  are  received  at  distributed  arrays.  In  this  paper,  we  investigate  new  approaches 
for  source  localization  with  multiple  arrays  that  exploit  partial  coherence  of  the  wavefronts  at 
distributed  arrays.  We  have  shown  previously  in  [1]  that  the  Cramer-Rao  lower  bound  (CRB)  on 
the  variance  of  source  location  estimates  is  reduced  when  coherence  from  array  to  array  is  exploited. 
We  also  showed  in  [1]  that  the  CRB  changes  little  for  suboptimum  source  localization  methods  that 
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employ  distributed  processing  to  reduce  the  communication  bandwidth  between  the  arrays  and  the 
fusion  center. 

In  this  paper,  we  focus  on  the  case  of  narrowband  processing,  and  algorithms  are  investigated 
for  source  localization  with  partially  coherent  signals  that  are  received  by  multiple  distributed 
sensor  arrays.  Two  algorithms  are  studied:  one  is  based  on  estimating  the  time  delay  in  the  signals 
received  by  widely-separated  arrays,  and  the  other  is  a  subspace  algorithm.  The  algorithms  are 
evaluated  in  terms  of  their  ability  to  achieve  the  performance  predicted  by  the  CRBs  developed  in 
[1],  and  by  processing  data  measured  at  Spesutie  Island,  Maryland.  Our  algorithms  and  testing 
address  the  issues  of  random  propagation  effects  that  result  in  partially  coherent  signals  arriving 
at  widely  separated  sensors,  and  Doppler  effects  due  to  source  motion  are  compensated  to  improve 
the  coherence. 

The  performance  analysis  presented  in  [1]  is  based  on  the  Cramer-Rao  lower  bound  (CRB)  on 
the  accuracy  of  any  unbiased  estimator  of  the  source  location.  We  model  the  signals  measured  at 
the  distributed  sensor  arrays  as  jointly  Gaussian  wideband  random  processes.  The  model  accounts 
for  propagation  effects  between  the  source  and  the  distributed  arrays,  including  frequency-selective 
spatial  coherence  and  different  signal  power  spectra  received  at  each  array.  The  spatial  coherence 
of  the  wavefronts  is  modeled  as  perfect  over  each  individual  array  but  imperfect  between  distinct 
arrays.  This  idealization  allows  us  to  study  the  effect  of  varying  coherence  between  arrays  on  source 
localization  accuracy.  Physical  modeling  of  frequency-selective  coherence  is  discussed  in  [15].  The 
power  spectral  density  of  the  source  is  arbitrary,  allowing  a  range  of  cases  to  be  modeled,  including 
narrowband  sources,  sums  of  harmonics,  and  wideband  sources  with  continuous  power  spectra. 

Previous  work  on  source  localization  with  acoustical  arrays  has  focused  on  angle  of  arrival  esti¬ 
mation  with  a  single  array  [2,  3,  4,  5].  These  works  use  the  coherent  wideband  focusing  approach 
[6,  7]  to  combine  correlation  matrices  from  different  narrowband  frequency  bins  into  a  single  corre¬ 
lation  matrix  that  admits  subspace  processing.  The  problem  of  imperfect  spatial  coherence  in  the 
context  of  narrowband  angle-of-arrival  estimation  with  a  single  array  has  been  studied  in  [8]-[ll]. 
Pauraj  and  Kailath  [8]  presented  a  MUSIC  algorithm  that  incorporates  the  nonideal  spatial  coher¬ 
ence,  assuming  that  the  coherence  variation  is  known.  Gershman  et  al.  [9]  provided  a  procedure  to 
jointly  estimate  the  spatial  coherence  loss  and  the  angles  of  arrival.  Song  and  Ritcey  [10]  provide 
maximum-likelihood  methods  for  estimating  the  parameters  of  a  coherence  model  and  the  angles 
of  arrival,  and  Wilson  [11]  incorporates  physical  models  for  the  spatial  coherence.  The  problem 
of  decentralized  array  processing  has  been  studied  in  [12]  and  [13].  Wax  and  Kailath  [12]  present 
subspace  algorithms  for  narrowband  signals  and  distributed  arrays,  assuming  perfect  spatial  coher¬ 
ence  across  each  array  but  neglecting  the  spatial  coherence  between  arrays.  Weinstein  [13]  presents 
performance  analysis  for  pairwise  processing  the  wideband  sensor  signals  from  a  single  array  and 
shows  negligible  loss  in  localization  accuracy  when  the  SNR  is  high. 

The  paper  is  organized  as  follows.  Section  2  describes  our  model  for  partially  coherent  signals 
observed  by  the  distributed  sensor  arrays.  Also  included  in  Section  2  are  results  from  measured 
data  to  illustrate  that  partial  signal  coherence  is  present  in  measured  aeroacoustic  data  with  sensors 
separated  by  hundreds  of  meters.  The  CRB  expressions  for  source  location  accuracy  are  reviewed 
in  Section  3.  Section  4  presents  fundamental  bounds  on  time-delay  estimation  that  are  an  extension 
of  the  Ziv-Zakai  bounds  [14]  to  the  case  of  partially  coherent  signals.  Section  5  presents  a  novel 
subspace  algorithm  for  narrowband  source  localization  with  distributed  arrays  and  partially  coher¬ 
ent  signals.  Examples  of  processing  measured  data  are  included  in  Sections  4  and  5,  and  Section  6 
contains  a  brief  discussion  of  our  continuing  work. 


Figure  1:  Geometry  of  source  location  and  H  distributed  sensor  arrays.  A  communication  link  is 
available  between  each  array  and  the  fusion  center. 


2  Data  Model 


A  model  is  formulated  in  this  section  for  the  discrete-time  signals  received  by  the  sensors  in 
distributed  arrays.  Consider  a  single  source  that  is  located  at  coordinates  (xs^yg)  in  the  (x,y) 
plane.  Then  H  arrays  are  distributed  in  the  same  plane,  as  illustrated  in  Figure  1.  Each  array 
h  £  H}  contains  Nh  sensors,  and  has  a  reference  sensor  located  at  coordinates  {xh,  yn)-  The 

location  of  sensor  n  G  {1, . .  • ,  Nh}  is  at  {xh  +  Axhn,  Vh  +  ^yhn)i  where  {Axhm  Ayhn)  is  the  relative 
location  with  respect  to  the  reference  sensor.  If  c  is  the  speed  of  propagation,  then  the  propagation 
time  from  the  source  to  the  reference  sensor  on  array  h  is 


-  Xh)^  +  (r/,  -  yh) 


(1) 


We  will  assume  that  the  wavefronts  are  well  approximated  by  plane  waves  over  the  aperture  of 
individual  arrays.  Then  the  propagation  time  from  the  source  to  sensor  n  on  array  h  will  be 
expressed  by  Th  +  Thm  where 
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--  [{cos<f>h)Axhn  +  (sin  (f>h)Ayhn] 


(2) 


where  Thn  is  the  propagation  time  from  the  reference  sensor  on  array  h  to  sensor  n  on  array  and 
(f>h  is  the  bearing  of  the  source  with  respect  to  array  h.  Note  that  while  the  far-field  approximation 
(2)  is  reasonable  over  individual  array  apertures,  the  wavefront  curvature  that  is  inherent  in  (1) 
must  be  retained  in  order  to  accurately  model  the  (possibly)  wide  separation  between  arrays. 

The  time  signal  received  at  sensor  n  on  array  h  due  to  the  source  will  be  represented  as  s/i(t  — 
Th  ~  Thn),  where  the  vector  of  signals  s{t)  =  [si  (t), . . . ,  Sif(t)]^  received  at  the  H  arrays  are  modeled 
as  real- valued,  continuous-time,  zero-mean,  wide-sense  stationary,  Gaussian  random  processes  with 
— oo  <  t  <  oo.  These  processes  are  fully  specified  by  the  H  X  H  cross-correlation  function  matrix 


Rsir)  =  E{s{t  +  T)s{tf}, 


(3) 


where  E  denotes  expectation,  superscript  T  denotes  transpose,  and  we  will  later  use  the  notation 
superscript  ♦  and  superscript  H  to  denote  complex  conjugate  and  conjugate  transpose,  respectively. 
The  (y,  h)  element  in  (3)  is  the  cross-correlation  function 


rs,gh{T)  =  E{sg{t  +  r)  Sh{t)} 


(4) 


between  the  signals  received  at  arrays  g  and  h.  The  correlation  functions  (3)  and  (4)  are  equivalently 
characterized  by  their  Fourier  transforms,  which  are  the  cross-spectral  density  functions 


/oo 

'^s,gh{T)  exp(-ia;r)  dr 
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(5) 


and  the  associated  cross-spectral  density  matrix 


G.(a;)  =  :F{R,(r)}. 


(6) 


The  diagonal  elements  Ga^hh{(*>)  of  (6)  are  the  power  spectral  density  (PSD)  functions  of  the  signals 
Sh{t),  and  hence  they  describe  the  distribution  of  average  signal  power  with  frequency.  The  model 
allows  the  average  signal  power  to  vary  from  one  array  to  another.  Indeed,  the  PSD  may  vary  from 
one  array  to  another  to  reflect  propagation  differences,  source  aspect  angle  differences,  and  other 
effects  that  lead  to  coherence  degradation  in  the  signals  at  distributed  arrays. 

Let  us  elaborate  the  deflnition  and  the  meaning  of  coherence  between  the  signals  Sg{t)  and 
Shit)  received  at  distinct  arrays  g  and  h.  In  general,  the  cross-spectral  density  function  (5)  can  be 
expressed  in  the  form 

(7) 

where  7a,jh(w)  is  the  spectral  coherence  function,  which  has  the  property  0  <  ^  1*  The 

coherence  function  'Ys^h{<^)  is  generally  complex-valued,  but  we  will  model  it  as  real-valued.  This 
is  a  reasonable  assumption  for  acoustic  propagation  environments  in  which  the  loss  of  coherence  is 
due  to  random  changes  in  the  propagation  path  length  [11, 15].  Note  that  our  assumption  of  perfect 
spatial  coherence  across  individual  arrays  implies  that  the  random  changes  in  the  propagation  path 
length  have  negligible  impact  on  the  intra-array  dela3re  in  (2).  Thus  the  apparent  bearing  <f>h  to  array 
h  remains  fixed,  but  the  inter-array  path  delays  are  randomly  perturbed  from  their  deterministic 
values  given  by  (1).  These  random  effects  are  modeled  by  the  signal  coherence  function,  'Ys,ghi<^)- 
We  model  the  signal  received  at  sensor  n  on  array  h  as  a.  sum  of  the  delayed  source  signal  and 
noise, 

Zhn (t)  =  Shit  -  Th  -  Thn)  +  Whn (0 ,  (8) 

where  the  noise  signals  Whnit)  are  modeled  as  real-valued,  continuous-time,  zero-mean,  wide-sense 
stationary,  Gaussian  random  processes  that  are  uncorrelated  at  distinct  sensors.  That  is,  the  noise 
correlation  properties  are 

E{v}gmit  +  =  »'«/(^)  ^ghSmn,  (9) 

where  r^uir)  is  the  noise  autocorrelation  function,  and  the  noise  power  spectral  density  is  Gt„(a;)  = 
We  then  collect  the  observations  at  each  array  h  into  x  1  vectors 
Zhit)  =  •  "}^h,Nkit)]^  for  h=  and  we  further  collect  the  observations  from  the  If 

arrays  into  a  (ATi  -1-  •  •  •  -f-  Nff)  x  1  vector 


The  elements  of  Z(t)  in  (10)  are  zero-mean,  wide-sense  stationary,  Gaussian  random  processes. 
We  can  express  the  cross-spectral  density  matrix  of  Z(t)  in  a  convenient  form  with  the  following 
definitions.  The  array  manifold  for  array  h  at  frequency  w  is 


a/i(w)  = 


expi-jtoThi) 

expi-jUThM) 


exp  [i^  ((cos  ^fc)  Aaifci  +  (sin  ^fc)  Ay^)] 
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(11) 


using  Thn  from  (2)  and  assuming  that  the  sensors  have  omnidirectional  response  to  sources  in  the 
plane  of  interest.  Let  us  define  the  relative  time  delay  of  the  signal  at  arrays  g  and  h  as 


^gh  —  ~ 


(12) 


where  is  defined  in  (1).  Then  the  cross-spectral  density  matrix  of  Z(t)  in  (10)  has  the  form 


Gz(w)  = 


(13) 


ai(w)ai(a;)^G,,ii(w) 
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Recall  that  the  source  cross-spectral  density  functions  Ga^gh{oj)  in  (13)  can  be  expressed  in  terms 
of  the  spectral  coherence  7a,jfc(w)  using  (7). 

Note  that  (13)  depends  on  the  source  location  parameters  (xa,ya)  through  a/i(a;)  and  Dgh. 
However,  (13)  points  out  that  the  observations  are  also  characterized  by  the  bearings  (f)i,. 
to  the  source  from  the  individual  arrays  and  the  relative  time  delays  Dgh  between  pairs  of  arrays. 
Therefore,  one  way  to  estimate  the  source  location  (®s,y»)  is  to  estimate  the  bearings  <^i, . 
and  the  pairwise  time  delays  Dgh. 


2.1  Signal  coherence  in  measured  data 

Next  we  present  results  from  measured  aeroacoustic  data  to  illustrate  typical  values  of  signal  co¬ 
herence  at  distributed  arrays.  The  experimental  setup  is  illustrated  in  Figure  2a,  which  shows  the 
path  of  a  moving  ground  vehicle  and  the  locations  of  four  microphone  arrays  (labeled  1,  3,  4,  5). 
Each  array  is  circular  with  N  =  7  sensors,  4-ft  radius,  and  six  sensors  equally  spaced  around  the 
perimeter  with  one  sensor  in  the  center.  We  focus  on  the  10  second  segment  indicated  by  the  O’s 
in  Figure  2a  (which  correspond  to  the  time  segment  340-350  sec  in  the  data).  Figure  2b  shows 
the  power  spectral  density  (PSD)  of  the  data  measured  at  arrays  1  and  3  during  the  10  second 
segment.  Note  the  dominant  harmonic  at  39  Hz.  Figure  2c  shows  the  estimated  coherence  between 
arrays  1  and  3  during  the  10  second  segment.  The  coherence  is  approximately  0.85  at  40  Hz,  which 
demonstrates  the  presence  of  significant  coherence  at  widely-separated  microphones.  Exploiting 
this  coherence  has  the  potential  for  improved  source  localization  accuracy.  Figure  2c  shows  the 
estimated  coherence  between  two  sensors  on  array  1,  spaced  by  8  feet.  Note  that  the  coherence  is 
close  to  unity  for  frequencies  in  the  range  from  about  40  to  200  Hz,  so  our  model  of  perfect  signal 
coherence  over  individual  arrays  seems  reasonable. 

The  Doppler  effect  due  to  source  motion  was  compensated  prior  to  the  coherence  estimate  shown 
in  Figure  2c.  Without  Doppler  compensation,  the  coherence  is  significantly  reduced,  as  shown  in 
Figure  3a.  The  time-varying  radial  velocity  of  the  source  with  respect  to  each  array  in  Figure  2a  is 
plotted  in  the  top  panel  of  Figure  3b.  If  s{t)  is  the  waveform  emitted  by  the  source  that  is  moving 
with  radial  velocity  v  with  respect  to  the  sensor,  then  the  sensor  receives  a  waveform  with  the  form 
s{at),  where  the  scaling  factor  a  is 

«  =  1  -  7  (14) 

and  c  is  the  speed  of  wave  propagation.  The  scaling  factor  a  is  plotted  in  the  bottom  panel  of 
Figure  3b.  Note  that  for  this  data  set,  0.98  <  a  <  1.02,  which  corresponds  to  a  Doppler  frequency 
shift  of  approximately  ±1  Hz  for  an  emitted  tone  at  50  hertz.  We  use  a  digital  resampling  algorithm 
to  compensate  for  the  Doppler  effect. 


Figure  2:  (a)  Path  of  ground  vehicle  and  array  locations  for  measured  data,  (b)  Mean  power 
spectral  density  (PSD)  at  arrays  1  and  3  estimated  from  measured  data  over  the  10  second  segment 
O  in  (b).  Top  panel  is  Gs^ii{f)^  bottom  panel  is  (c)  Mean  spectral  coherence  js.isif) 

between  arrays  1  and  3  estimated  over  the  10  second  segment,  (d)  Mean  spectral  coherence  for  two 
sensors  on  array  1,  with  sensor  spacing  8  feet. 
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Figure  3:  (a)  Mean  short-time  spectral  coherence  between  arrays  1  and  3  if  Doppler  is  not  compen¬ 
sated.  (b)  Radial  velocity  and  Doppler  scaling  factor  a  in  (14)  for  source  in  Figure  2a  with  respect 
to  each  array. 


3  CRBs  on  Localization  Accuracy 


The  problem  of  interest  is  to  estimate  the  source  location  parameter  vector  ©  =  [xs,ys]^  using  T 
samples  of  the  sensor  signals  Z(0),  Z{Ta), . .  ■,Z{{T  —  1)  -Tg),  where  Tj  is  the  sampling  period.  Let 
us  denote  the  sampling  rate  by  /,  =  l/T^  and  Ua  =  27r/s.  We  will  assume  that  the  continuous¬ 
time  random  processes  Z{t)  are  band-limited,  and  that  the  sampling  rate  fa  is  greater  than  twice 
the  bandwidth  of  the  processes.  Then  Priedlander  [18]  has  shown,  using  a  theorem  of  Whittle 
[19],  that  the  Fisher  information  matrix  (FIM)  J  for  the  parameters  0  based  on  the  samples 
Z(0),  Z{Ta), . . Z((T  —  1)  •  Ta)  has  elements 


J-  =  ——  f 

2Ua  Jo 


\  d9i 


Gziu))' 


89, 


GzM  i,j  =  1,2, 


(15) 


where  “tr”  denotes  the  trace  of  the  matrix.  The  CRB  matrix  C  =  J“^  then  has  the  property 
that  the  covariance  matrix  of  any  unbiased  estimator  0  satisfies  Cov(©)  —  C  >  0,  where  >  0 
means  that  Cov(©)  —  C  is  positive  semidefinite  [17].  The  CRB  provides  a  lower  bound  on  the 
performance  of  any  unbiatsed  estimator.  Equation  (15)  provides  a  convenient  way  to  compute  the 
FIM  for  the  distributed  sensor  array  model.  It  provides  a  powerful  tool  for  evaluating  the  impact 
that  various  parameters  have  on  source  localization  accuracy.  Parameters  of  interest  include  the 
spectral  coherence  between  distributed  arrays,  the  signal  bandwidth  and  power  spectrum,  the  array 
placement  geometry,  and  the  SNR.  The  FIM  in  (15)  is  not  easily  evaluated  analytically,  but  it  is 
readily  evaluated  numerically  for  cases  of  interest. 

Consider  an  acoustic  source  that  has  a  narrowband  power  spectrum.  That  is,  the  PSD  Gs^hh{^) 
of  the  signal  at  each  array  h  =  1, . .  is  nonzero  only  in  a  narrow  band  of  frequencies  ujq  — 


(Aw/2)  <  w  <  Wo  +  (Aw/2).  If  the  bandwidth  Aw  is  chosen  small  enough  so  that  the  w— dependent 
quantities  in  (15)  are  well  approximated  by  their  value  at  wq,  then  the  narrowband  approximation 
to  the  FIM  (15)  is  ’ 


TAw^  /aGz(wo) 
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dOi 


Gz(wo) 
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The  quantity  2^  multiplying  the  FIM  in  (16)  is  the  time-bandwidth  product  of  the  observations. 
In  narrowband  array  processing,  the  T  time  samples  per  sensor  are  often  segmented  into  M  blocks 
containing  T/M  sample  each.  Then  the  discrete  Fourier  transform  (DFT)  is  applied  to  each  block, 
and  the  complex  coefficients  at  frequency  cjq  (at  each  sensor)  are  used  to  form  M  array  “snapshots” . 
In  this  case,  the  quantity  is  approximately  equal  to  M. 


4  Threshold  Coherence  for  Time  Delay  Estimation 

In  this  section,  we  present  a  model  for  partial  signal  coherence  in  terms  of  an  equivalent  additive 
noise  component,  which  allows  the  fundamental  bounds  on  time  delay  estimation  in  [14]  to  be 
extended  to  the  case  of  partially  coherent  signals.  Bounds  on  time  delay  estimation  are  important 
because  one  of  the  suboptimum  processing  methods  proposed  in  [1]  for  distributed  processing  with 
low  communication  bandwidth  requires  time  delay  estimation  between  widely  separated  sensors. 
The  fundamental  bounds  in  [14]  are  useful  because  they  specify  the  required  signal-to-noise  ratio 
(SNR)  such  that  the  CRB  is  attainable.  That  is,  if  the  SNR  is  less  than  a  threshold,  then  the 
CRB  on  source  localization  accuracy  is  an  optimistic  and  unattainable  lower  bound.  We  will  show 
that  for  signals  that  are  partially  coherent  when  received  at  the  sensors,  a  threshold  phenomenon 
occurs  with  respect  to  coherence.  That  is,  the  signal  coherence  must  exceed  a  threshold  in  order 
for  the  CRBs  to  be  attainable.  The  existence  of  a  threshold  phenomenon  for  signal  coherence  is  a 
useful  refinement  of  [1],  since  the  CRBs  in  [1]  are  formulated  without  regard  for  the  conditions  of 
attainability. 

The  appendix  presents  a  decomposition  of  partially  coherent  signals  into  a  coherent  component 
and  an  incoherent  additive  noise  component.  Consider  estimation  of  the  time  delay  D  in  the  model 
presented  as  (54), (55)  and  (64), (65)  in  the  appendix.  Let  us  specialize  to  narrowband  processing, 
with  signal  bandwidth  Aa;  centered  at  ujo  and  observation  time  T  seconds.  Further,  we  assume  that 
the  signal  power  is  identical  at  each  sensor,  and  we  define  the  following  constants  for  notational 
simplicity: 

Gs, =  G^5,2(^o)  =  Gs^  Gw{(^o)  =  Gwy  Ts, 12(^0)  =  Js*  (17) 

Then  combining  the  formulation  in  the  appendix  with  the  development  in  [14],  we  can  show  that 
the  following  SNR  expression  characterizes  the  performance  of  time  delay  estimation  with  partially 
coherent  signals: 
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From  [14],  the  threshold  SNR  for  CRB  attainability  in  the  narrowband  time  delay  estimation 
problem  is 


SNRthresh  — 


24a;2  ) 


(19) 


(20) 


where  <f>{y)  =  Jy/^  exp(-t^/2)  dt.  Thus 

SNR(7,)  >  SNRthresh 


identifies  the  values  of  signal  coherence  jg  and  signal/noise  PSDs  GgfGw  for  CRB  attainability. 
We  can  combine  (18)  with  (20)  to  obtain  the  condition  for  CRB  attainability 


1 


(21) 


which  is  possible  only  if 

^  r-j — —1 — • 

^  f  SNRthr«h 

For  a  specific  narrowband  time  delay  estimation  scenario,  the  threshold  SNR  for  CRB  attainability 
is  given  by  (19),  and  (22)  provides  a  corresponding  threshold  coherence  for  CRB  attainability. 

For  example,  consider  time  delay  estimation  in  a  band  centered  at  ljq  =  27r50  rad/sec  using  a 
time  segment  of  duration  2  seconds. 


•  For  a  1  Hz  bandwidth,  i.e.,  Au;  =  2n  rad/sec,  nearly  perfect  coherence  (7,  =  0.99995)  is 
required  in  order  to  enter  the  regime  of  CRB  attainability. 


For  a  10  Hz  bandwidth,  i.e  Au  =  27rl0  rad/sec,  the  coherence  must  exceed  0.93  for  CRB 
attainability. 


The  implication  is  that  for  the  narrowband  signals  with  Au?  =  2it  rad/sec,  any  loss  of  coherence 
leads  to  difficulty  with  time  delay  estimation  due  to  ambiguities  arising  from  the  narrowband  signals. 
But  for  sufficiently  wideband  signals,  e.g.,  Au  =  27rl0  rad/sec,  a  certain  amout  of  coherence  loss 
can  be  tolerated  while  still  allowing  unambiguous  time  delay  estimation. 

We  present  an  illustration  based  on  processing  the  measured  data  for  the  source  in  Figure  2a. 
Figure  4  shows  results  of  cross-correlation  processing  of  the  data  for  a  2  second  segment  at  time 
342  seconds.  Figure  4a  is  obtained  by  cross-correlating  the  signals  received  at  arrays  1  and  3,  for 
which  the  coherence  is  appreciable  only  over  a  narrow  band  near  39  Hz.  A  peak  in  the  cross¬ 
correlation  is  not  evident,  which  is  expected  based  on  the  preceding  analysis,  since  nearly  perfect 
coherence  is  needed  for  narrowband  time  delay  estimation  in  this  scenario.  Figure  4b  is  obtained 
by  cross-correlating  the  signals  received  at  two  sensors  on  array  1,  where  the  coherence  is  as  shown 
in  Figure  2d.  Due  to  the  high  signal  coherence  over  a  wide  band,  a  peak  is  clearly  evident  in  the 
cross-correlation . 


5  Subspace  Processing 

In  this  section,  we  begin  with  an  eigenanalysis  of  the  cross-spectral  density  matrix  (13)  for  the 
case  of  =  2  arrays  containing  Ni  =  N2  —  N  sensors  operating  in  a  narrow  band  of  frequencies 
centered  at  uq.  This  analysis  leads  to  a  MUSIC-like  subspace  algorithm  for  source  localization 
with  distributed  arrays  and  partial  signal  coherence.  We  illustrate  with  computer  simulations 
and  measured  aeroacoustic  data  that  the  new  subspace  algorithm  is  limited  in  performance  due  to 
source  location  ambiguities  that  arise  from  the  large  separation  between  arrays  and  the  narrowband 
signals. 

To  simplify  notation,  we  will  let  SLh  represent  the  array  manifold  a;i(a;o),  represent  the  average 
signal  power  at  array  h  in  the  frequency  band  of  interest,  is  the  average  noise  power,  and  7  is 
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CROSS-CORRELATION:  ARRA  VI ,  SB^SORS 


Figure  4:  (a)  Cross-correlation  of  signals  at  arrays  1  and  3  in  Figure  2a  for  source  at  time  342  sec 
(b)  Cross-correlation  of  signals  at  two  sensors  on  array  1,  separated  by  8  feet. 


the  coherence  12(^0)-  We  allow  7  to  be  complex- valued  with  I7I  <  1.  Then  (13)  can  be  expressed 
as  the  correlation  matrix 


R(12)  = 


TO’iO’2 


T  0-1 0-2 


1 


<S)  InxN  ° 


aiaf 

a2  af  exp(+ja>oI>i2) 


ai  a^  exp{-ju)oDu) 
a2af 


+ 


(23) 

where  ®  denotes  Kronecker  product,  o  denotes  elementwise  product,  and  livxiv  is  an  iV  x  TV  matrix 
of  I’s.  The  superscript  (12)  in  (23)  indicates  that  the  quantity  is  the  combined  correlation  matrix 
of  all  sensors  in  arrays  1  and  2.  The  correlation  matrices  for  the  individual  arrays  1  and  2  are 


(rlakB-h  +  h=l,2. 


(24) 


It  is  well-known  that  the  individual  array  correlation  matrices  and  in  (24)  have  a  rank-1 
signal  subspace,  and  that  the  signal  subspace  eigenvalue  is 

\^^^  =  Nal  +  cTl,  h=l,2 

with  corresponding  eigenvector 

eW=  a;,,  h=  1,2.  (26) 

The  remaining  eigenvalues  n  =  2, . . N  are  equal  to  cr^. 

Now  we  present  the  eigenanalysis  of  the  combined  correlation  matrix  R(^^)  in  (23).  The  deriva¬ 
tion  of  these  results  is  similar  to  the  reasoning  in  the  appendix  of  [8].  We  assume  that  the  arrays 
are  labeled  so  that  >  <^1  and  that  the  sensors  are  omnidirectional  so  that  the  array  manifold 


elements  are  of  the  form  [a^Jn  =  exp(-ywoTXft).  The  rank  of  the  signal  subspace  component  of 
is  at  most  two,  and  the  two  largest  eigenvalues  are 


The  corresponding  eigenvectors  are 


(27) 
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(30) 


provided  Some  observations  are  as  follows. 

•  In  the  case  that  <ti  =  <T2  =  <r,,  (27)-(30)  reduce  to 

A<“>  =  W<^(l  +  |Tl)+ffi,  4“>  =  N,^(l-|7|)  +  oii, 


(31) 

(32) 


Note  that  (32)  is  vaUd  as  long  as  ^  ^  0,  and  in  this  case  the  dgenvectors  are  independent  of 
the  coherence  ma^tnde.  In  the  case  of  high  SNR  so  that  «<  “  ItDi  tli®  of 

signal  dgenvalues  is  approximately 
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f^om  (33),  in  order  for  so  that  A[^*^  is  the  dominant  mgenvalue,  the  coherence 

must  satisfy  fy]  >  0.8181. 

•  In  the  case  of  perfect  coherence  fy(  =  1,  R(i^)  has  a  rank-1  signal  subspace,  with  signal 
eigenvalue  and  mgenvector 


^{12, OOH)  ^  ^  ^2)  ^  ^2  ^  e(12.COH)  ^ 
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(34) 


The  expressions  in  (34)  are  the  natural  extension  of  (25)  and  (26)  applied  to  the  *^uperarray” 
of  all  sensors,  allowing  for  different  signal  power  levels  at  the  two  arrays. 

•  In  the  case  of  partial  coherence  H  <  1,  R(^^)  has  a  rank-2  signal  subspace.  The  signal 
dgenvalues  are  ordered  A[^*^  >  A^'*  ,  with  equality  if  and  only  if  7  =  0  and  oj  =  02*  In  the 
incoherent  case  7  ==  0,  the  mgnal  dgenvectors  are 
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(35) 


Thus  in  the  incoherent  case,  the  signal  eigenvectors  (35)  of  are  completely  determined 
by  the  signal  eigenvectors  (26)  of  individual  array  correlation  matrices  and  R^^^. 


Note  that  the  span  of  the  R^*^^  signal  subspace  eigenvectors  in  (29)  and  (30)  is  equal  to  the 
span  of  the  incoherent  eigenvectors  in  (35).  Therefore,  for  I7I  <  1,  the  signal  subspace  of  R^^^) 
is  spanned  by  the  signal  subspace  eigenvectors  of  the  individual  array  correlation  matrices  R^^^ 

r  0 

and  It  follows  that  R(^^)  and  the  block  diagonal  matrix 


0  R(2)  ‘  identical 

signal  and  noise  subspaces.  Thus  a  subspace  algorithm  that  exploits  joint  processing  of  the 

in 


.(12) 


two  arrays  must  exploit  the  structure  of  the  dominant  signal  subspaee  eigenvector  e 
(29). 

The  following  subspace  procedure  may  be  used  to  estimate  the  source  location  (xsjJ/s)  using 
the  correlation  matrices  R(^1,  RI^I,  and 

1.  Form  the  correlation  matrices  r1^)  and  R^^)  using  the  observations  from  the  individual  arrays. 
Elstimate  the  noise  power  <r^  from  the  noise  eigenvalues  of  R^^^  and  R(^^.  Then  estimate  the 
signal  power  at  each  array  using  the  largest  eigenvalue  of  R(^)  and  R^^^  with  (25) : 


=  h=l,2. 


(36) 


2. 


Form  the  joint  correlation  matrix  and  estimate  the  coherence  magnitude  I7I  using  the 


two  largest  eigenvalues  A^^^^  of  R^^^)  with  (27)  and  (28): 


(37) 


as  long  as  <ri  ^  <r2*  Note  that  if  <ri  =  <r2j  then  it  is  not  necessary  to  estimate  the  coherence 
magnitude  (see  (32)  and  step  3  below).  We  vdll  assume  that  7  is  real-valued  and  positive, 
which  is  motivated  by  physical  considerations  as  discussed  earlier  in  this  paper  and  in  [11, 15]. 


3. 


Use  the  estimates  for  ai,<7-2,  and  I7I  in  the  signal  eigenvector  expression  for  in  (29)  to 
create  a  “search  vector”  g: 


g(*s.y*) 
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2  \<r2  <ri/ 

ai 

T* 
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(38) 


where  ai,a2,  and  D12  are  functions  of  the  source  location  (®»,y«).  (If  ci  =  0'2,  then  the 
search  vector  is  based  on  in  (32)  instead.) 


4.  Let  Cn  n  =  1, . . . ,  2N  denote  the  eigenvectors  of  ordered  such  that  the  corresponding 
eigenvalues  follow  A^^^^  >  •  •  •  >  Then  form  the  matrix 

v=[4‘" 

A  MUSIC-type  spectrum  may  then  be  defined  as 

Pc{xt,  y,)  =  gHvv^g  ’ 

and  the  value  of  (x,,  y,)  that  maximizes  (40)  is  an  estimate  of  the  source  location. 


5.  Wax  and  Kailath  [12]  proposed  an  “incoherent”  MUSIC-type  source  location  estimator  that 
ignores  coherence  between  the  arrays.  The  form  of  the  Wax  and  Kailath  estimator  can  be 
expressed  as  follows, 


where  •  •  •  ?  contains  the  noise  eigenvectors  of 

This  analysis  extends  to  the  case  of  H  arrays  with  a  single  source.  If  the  signals  are  partially 
coherent  at  all  arrays,  then  the  dimension  of  the  signal  subspace  of  the  combined  correlation  matrix 
is  H.  If  K  signals  are  present,  then  the  combined  correlation  matrix  has  signal  subspace  with  rank 
KH. 

5.1  Examples  of  subspace  processing 

We  present  examples  of  the  “partially  coherent”  MUSIC  algorithm  (40)  and  the  “incoherent” 
MUSIC  algorithm  (41)  with  computer  simulations  and  measured  data.  First,  we  simulate  a  scenario 
with  H  =  2  arrays.  The  individual  arrays  are  identical  and  contain  =  JV2  =  TV  =  7  sensors. 
Each  array  is  circular  and  has  4-ft  radius,  with  six  sensors  equally  spaced  around  the  perimeter  and 
one  sensor  in  the  center.  Narrowband  processing  centered  at  40  Hz  is  assumed,  with  an  SNR  of  40 
dB  per  sensor  at  array  1  and  34  dB  at  array  2,  i.e.,  Gs^ii{u)/Gw{u)  =  10^,  (7«,22(^)/Gty(^)  =  2500. 
The  coherence  is  set  to  two  values,  7,^12  =  0.9  and  0.999.  The  arrays  are  located  at  coordinates 
(xi,yi)  =  (0,150),  (x2,i/2)  =  (0,0),  and  one  source  is  located  at  (x^ji/^)  =  (200,300),  where  the 
units  are  meters.  Source  location  estimation  is  performed  using  M  =  30  snapshots,  where  each 
snapshots  contains  an  estimate  of  the  measured  complex  amplitude  at  40  Hz  at  each  sensor. 

Results  are  presented  in  Figure  5,  where  part  a  shows  an  approximate  CRB  as  a  function  of 
signal  coherence.^  Increasing  coherence  leads  to  modest  improvement  in  the  localization  potential. 
Figures  5b  and  5c  contain  representative  spectra  of  incoherent  and  partially  coherent  MUSIC, 
respectively,  for  coherence  7,^12  =  0.9.  Note  that  the  incoherent  MUSIC  spectrum  has  a  rather 
broad  but  unique  maximum.  The  partially  coherent  MUSIC  spectrum  has  a  sharp  ridge  in  one 
direction,  with  multiple  peaks  along  the  ridge.  The  multiple  peaks  arise  due  to  the  ambiguities  in 
source  location  caused  by  the  narrowband  signals  and  large  separation  between  arrays. 

Table  1  contains  results  of  Monte  Carlo  simulations  with  100  runs  for  both  coherence  values, 
7«,12  =  0.9  and  0.999.  The  incoherent  MUSIC  performs  close  to  the  approximate  CRB  in  Figure  5a 
corresponding  to  75^12  =  0.  Partially  coherent  MUSIC  performs  worse  than  incoherent  MUSIC 
because  of  the  ambiguities:  partially  coherent  MUSIC  sometimes  chooses  the  incorrect  peak  along 
the  ridge  in  the  spectrum  shown  in  Figure  5c.  Note  that  partially  coherent  MUSIC  is  more  accurate 
for  7,  =  0.999  than  for  7,  =  0.9,  suggesting  that  the  ambiguities  are  less  severe  when  the  signal 
coherence  is  larger. 

We  have  applied  incoherent  MUSIC  and  partially  coherent  MUSIC  to  the  measured  data  sce¬ 
nario  in  Figure  2.  The  data  is  processed  in  a  narrow  band  around  39  Hz  for  the  2-second  interval 
between  345  and  347  seconds.  Both  algorithms  produce  identical  estimates  for  the  source  location 
in  this  case,  but  the  spectra  for  partially  coherent  MUSIC  are  shown  in  Figure  6.  The  top  panel 
clearly  shows  the  ambiguities  that  limit  the  performance  with  narrowband  signals.  The  contour 
plot  in  the  bottom  panel  also  shows  the  periodic  structure  of  the  ambiguities. 

'The  CRB  is  approximate  because  the  mapping  of  the  scaling  factor  in  (16)  to  the  number  of  snapshots  M 
is  approximate. 
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Figure  5:  (a)  CRBs  on  localization  accuracy,  (b)  Incoherent  MUSIC  and  (c)  partially  coherent 
MUSIC  spectra. 


Coherence 

Approx.  CRB 

Incoherent  MUSIC 

Partially  Coherent  MUSIC 

Xs 

Vs 

Xs 

Vs 

X.n 

Vs 

0.9 

1.17 

1.10 

1.36 

1.10 

2.14 

1.71 

0.999 

0.81 

0.87 

1.41 

1.16 

1.96 

1.61 

0.0 

1.50 

1.22 

Table  1;  RMS  error  in  simulation  results  for  incoherent  and  partially  coherent  MUSIC. 
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Figure  6:  Partially  coherent  MUSIC  spectrum  for  measured  data. 


6  Concluding  Remarks 

Items  of  continuing  work  include  the  following, 

•  The  analysis  in  this  paper  was  focused  primarily  on  narrowband  processing.  We  found  that 
source  location  ambiguities  due  to  the  narrowband  signals  are  an  important  limitation  when 
distributed  arrays  are  processed  jointly.  Signals  that  maintain  coherence  over  a  wider  band¬ 
width  are  expected  to  experience  fewer  problems  with  ambiguities,  so  we  are  evaluating 
algorithms  for  this  scenario. 

•  We  are  continuing  to  process  measured  data,  particularly  data  from  distributed  arrays  with 
synchronized  sampling  at  the  arrays,  courtesy  of  Sanders  Corporation. 

•  We  are  investigating  models  for  the  source  motion  that  will  improve  accuracy  and  allow 
tracking. 

•  We  are  investigating  extensions  to  the  case  of  tracking  multiple  moving  sources. 


Appendix 


A  Equivalent  Additive  Noise  Model  for  Coherence 

The  formulation  in  this  appendix  begins  with  Gaussian  random  variables  for  simplicity,  and  then 
the  formulation  is  extended  to  Gaussian  random  processes. 

Let  X  and  Y  be  two  complex,  circular,  Gaussian  random  variables  with  zero  mean  and 

E{\X\^}  =  <rx,E{\Y\^}  =  4.  E{XY*}  =  pcxCY  (42) 


Then  X  and  Y  can  be  expressed  as 


X 

Y 


axS  +  Ni 


l/>l 


,-oy5  +  Ni, 


(43) 

(44) 


where  S,  N\,  Ni  are  independent,  zero  mean,  complex,  circular  Gaussian  random  variables  with 

E{m  =  \p\  (45) 

=  4(i-H)  (46) 

E{\Ni\^}  =  (47) 


The  relations  in  (42)  can  be  verified  using  (43)- (47).  Note  that  for  X,  the  “coherent”  part  crxS 
has  variance  1/>|4)  while  the  “incoherent”  part  iVi  has  variance  (1  —  \p\)o\‘  As  |p|  decreases,  a 
larger  portion  of  4  applied  to  N-y ,  so  the  incoherent  noise  gets  stronger. 

Suppose  that  X  and  Y  described  by  (42)  are  observed  in  noise. 


Zy  =  X  +  Wy  =  (rxS  +  Ny  +  Wy 

Zi  =  Y  +  Wy  =  fr<TYS  +  Ni+Wi, 
\P\ 


(48) 

(49) 


where  Wi  and  W2  are  independent,  zero  mean,  complex,  circular  Gaussian  random  variables  with 
variance  cr^,  and  Wi,W2  are  independent  of  X,  Y  (and  therefore  S,  Ni,  N2).  If  we  define  the  SNR 
of  Zi  and  Z2  based  only  on  the  coherent  signal  component,  then 


SNRi 

SNR2 
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(50) 

(51) 


Some  observations  about  (50)  and  (51)  follow. 

1.  SNRi  and  SNR2  increase  monotonically  as  |/9|  increases. 


2.  SNRi  (SNR2)  increases  monotonically  as  (cry)  increases,  and  for  fixed  \p\  is  bounded  by 

lim  SNRi=  lim  SNR2  =  J47.  (52) 

<7^-fOO  1  “  \p\ 


Thus  a  loud  source  (large  (Jx^^y)  larger  coherent  SNR,  i.e.,  increasing  cannot 

reduce  the  coherent  SNR,  but  the  limit  (52)  cannot  be  exceeded  for  a  given  \p\. 


3.  If  a  particular  coherent  SNR  is  desired,  then  (52)  implies  a  threshold  coherence  magnitude 
\p\  that  is  required  to  achieve  that  SNR,  given  by 


,  SNR 
-  l  +  SNR' 


(53) 


Note  that  for  large  coherent  SNR,  nearly  perfect  coherence  \p\  «  1  is  required. 

The  analysis  extends  to  complex  Gaussian  random  processes  as  follows.  Consider  a  time  delay 
estimation  problem  with  two  sensors  and  noisy  observations  of  the  form  in  (8), 

Zi{t)  =  Si{t)+Wi{t)  (54) 

Z2{t)  =  S2{t- D)  +  W2{t),  (55) 

where  t«i(t),  u;2(t)  are  additive,  white,  Gaussian  noise  (AWGN)  processes,  si(t),S2(^)  are  partially 
coherent  Gaussian  random  processes,  and  D  is  the  time  delay.  The  AWGN  processes  wi{t),W2{t) 
have  properties  as  in  (9)  with  power  spectral  density  (PSD)  Gt„(w),  and  the  signals  si(t),S2(0  are 
characterized  by  the  cross-spectral  density  matrix 

t.  .\  tr^  I.  .\r^  I.  .\\l/2  1 

(56) 


(57) 

(58) 


Sl(t) 

Ga,ii(a;) 

Tj>,12(‘*’)  {Gs,\l{<^)Ga,22{<^))^^^ 

S2(t) 

.  T»,12(w)*  (Ga,ii(w)Ga,22(w))^^^ 

Gs,22{^) 

Then  analogous  to  (43,  (44),  si(f)  and  S2{t)  can  be  represented  as 

Si  (t)  =  hi  {t)  *s{t)  +  ni  (t) 

S2(t)  =  /i2(0  *s{t)+  n2{t) 


where  s{t),ni{t),n2{t)  are  independent,  zero  mean,  circular,  complex  Gaussian  random  processes, 
♦  denotes  convolution,  and 


Hi(pj)  =  (59) 

=  r^^4^G..22(w)'/'  (60) 

It.,12(w)| 

G,(w)  =  |7,,i2(w)1  (61) 

Gi(a>)  =  G,,ii(a;)  [1  -  |7,,i2(w)0  (62) 

^2(0;)  =  G*,22(w)  [1  —  1t»,12(‘*>)|]  •  (63) 


Ga(tj),  (j2(w)  are  the  PSDs  of  s(t),  ni(t),  n2(t),  and  H\  {u),H2{(jj)  are  the  frequency  responses 

corresponding  to  hi{t),h2{t)  that  model  the  deterministic  “channels”  from  source  to  sensors.  Then 
(57), (58)  can  be  inserted  into  (54), (55), 

zi  {t)  =  {hi  *  s)  (t)  +  ni  (t)  +  wi  {t)  (64) 

^2(0  =  (^2  *  s){t  -  D)  +  n2{t)  +  tU2(0)  (65) 

which  fits  the  standard  model  for  time  delay  estimation  of  coherent  signals  observed  through  linear 
filters  in  AWGN  [20].  The  partial  signal  coherence  7s,i2(w)  between  si{t),S2{t)  in  (54), (55)  is 
equivalently  modeled  by  the  filtered  coherent  signals  {hi  *  s) ,  (^2  *  s)  and  the  excess  additive  noise 
ni{t),n2{t). 
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